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A method  for  direct  integration  of  the  dynamic  governing  partial  differential 
equations  of  motion  for  structural  members  is  presented.  This  technique  is  • 
called  the  continuous  space  discrete  time  (CSDT)  Riccati  transfer  matrix 
method.  Numerical  results  for  bar  and  beam  example  problems  indicate  that  the 
method  is  numerically  stable  and  accurate  for  calculating  the  dynamic  response 
of  linear  structural  members. 
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INTRODUCTION 


The  most  common  approach  to  structural  dynamics  problems  is  to 
discretize  the  continuous  structure  in  space,  using,  for  example,  the 
finite  element  model.  This  leads  to  a set  of  ordinary  differential  equa- 
tions in  time.  Integration  schemes  or  modal  superposition  are  employed 
to  solve  these  equations. 

In  the  present  effort  the  governing  partial  differential  equations 
for  structural  members  such  as  rods,  plates,  and  shells  are  transformed 
into  ordinary  differential  equations  in  space  by  discretizing  the  time 
derivatives  using  a finite  difference  scheme.  A mixed  method  for  struc- 
tural members,  such  as  the  transfer  matrix  method  for  beams  or  numerical 
integrations  for  shells,  is  then  applied  to  solve  these  ordinary  diff- 
erential equations.  These  computations  in  space  can  be  stabilized  with  the  aid 
of  the  field  method  or  Riccati  transformations. 

The  method  developed  in  the  present  paper  is  called  the  continuous- 
space  discrete  - time  (CSDT)  Riccati  transfer  matrix  method  since  only 
the  time  variable  is  treated  in  discrete  form  and  the  Riccati  transfer 
matrix  method  Cl3  is  employed  to  eliminate  the  numerical  instabilities 
so  often  encountered  C 2]  in  spatial  calculations.  As  mentioned  above, 
this  method  differs  from  the  commonly  used  direct  integration  method  for 
structural  dynamics,  in  which  the  numerical  Integration  is  performed  on 
a system  of  second  order  differential  equations  resulting  from  the  usual 
structural  approximations  of  the  spatial  geometry  of  the  members  C3j. 

The  continuous- space  discrete- time  method  has  been  used  in  analog 
or  hybrid  computers  for  the  solution  of  partial  differential  equations 
C^»  5,  6,  7,  8 ].  The  central  difference  formulation  was  used  in  discretizing 
the  time  derivatives.  Stability  of  the  method  in  the  direction  of  time 
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has  been  studied  for  parabolic  heat  transfer  equations  [9]  and  two 
methods,  i.e.  the  method  of  decomposition  [10]  and  the  invariant 
imbedding  method  [11,12],  were  developed  to  eliminate  the  instability 
associated  with  the  resulting  boundary  value  problem  for  the  ordinary 
differential  equations  in  space.  Breed  [13]  used  a form  of  the  CSDT 
approach  for  the  transient  analysis  of  rotating  shafts.  A central 
difference  time  discretization  was  used  in  his  formulation.  However, 
no  effort  was  made  to  stablize  the  integration  for  the  resulting  spatial 
boundary  value  problem. 

Without  such  a stabilization,  the  applications  are  severely  limited. 
This  paper  extends  the  continuous-space  discrete-time  method  to  the 
transient  analysis  of  structural  members.  The  Newmark  generalized 
acceleration  formulation  is  used  for  time  discretization.  Advantages 
of  this  formulation  over  the  central  difference  formulation  usually 
adopted  for  the  continuous-space  discrete-time  method  will  be 
demonstrated. 


ACCESSION  tm  . 

NT IS  Wkiw  Sactfcn 

DOC  Bid*  lection 

UtWNNP-TWB 


JVSn.TCAT 

•Tif 

DISUira/HYXillSBIinTf  (MS 

Disl.  A'Ail  j.’d/ur  SPECIAL 

ft 

II.  3 


□ □ 


FORMULATION  OF  THE  CONTINUOUS- SPACE  DISCRETE-TIME 
RICCATI  TRANSFER  MATRIX  METHOD 


In  general,  the  partial  differential  equations  of  motion  for 
structural  members  can  be  expressed  as 


DW(x, t) 


j h Aj(x) 


3jW(x,t) 

3tj 


-F(x,t) 


(1) 


where  W(x,t)  is  the  N-dimensional  column  vector  of  dependent  variables, 

Aj (x)  is  an  N-square  spatial  matrix,  D(x)  is  an  N-square  spatial  matrix 
linear  differential  operator  and  F(x,t)  is  the  N-dimensional  force  vector. 
The  terms  associated  with  A^  and  A^  are  usually  identified  as  the  damping 
and  inertia  terms  respectively.  If  some  form  of  finite  difference 
discretization  is  used  to  approximate  the  inertia  and  damping  terms  in 
equation  (1),  then  the  partial  differential  equations  can  be  transformed 
into  an  ordinary  differential  equation  with  spatial  derivatives  only. 

The  linear  multistep  discretization  is  used  for  this  purpose  due  to  the 
fact  that  most  of  the  commonly  used  direct  integration  methods  in 
structural  dynamics,  for  example,  the  central  difference  method  ClO, 
the  Houbolt  method  C^3,  the  Newmark  method  and  the  stiffly  stable 

method  C 173,  can  be  derived  from  the  linear  multistep  formulation.  This 
formulation  can  be  written  as 


1fQ  “i  WiH-I-i  " At  i-0  6i  ^n+l-i 

(2) 

P 9 P .. 

Z y W > At  I 6 W 

i-0  ' i n+l-i  i-0  °i  n+l-i 

(3) 

where  dot  notation  is  used  to  represent  time  derivatives,  ct^,  0^,  y^  and 
6i  are  coefficients,  At  is  the  time  step  and  the  subscript  n+l-i  denotes 

the  time  at  (n+l-i)At  with  n-0,  1,  . Rearrange  equations  (2)  and  (3) 

in  the  form 


II. 4 


VJ  S — TJ  + P 

n+1  B At  n+1  rn 
o 

w - — — w + ^ 

n+1  6 At2  n+1 
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P 

where  P - — 7 a , W , . . 

n S it  /.  i n+l-i  8 
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6 At2  i*l 

o 
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The  functions  Pn  and  Qn  involve  variables  at  previous  times  only  and 
hence  can  be  considered  as  the  historical  part  in  the  formulation. 
Substituting  equations  (4)  and  (5)  into  the  governing  equations  of 
motion  (1)  at  time  (n+1) At,  we  have 

DVl  ' (Fit  + TTT'Vl  + Vn  * *A  - Fn+1  <8) 

o 6 At^ 
o 

°r  DWn+l  - ^n+l  " Rn  <9> 

A a Ay 

where  K - — — + 7=777  (10) 


Rn  * Vn  + A2^n  ' Fn+1  (11) 

Equation  (9)  is  a differential  equation  with  spatial  derivatives  only 

and  hence,  the  dynamic  analysis  of  the  structural  members  can  be  treated 

at  each  time  step  (n-0,1, ) as  a static  problem  by  considering  R to 

n 

be  a generalized  external  force  acting  on  the  member. 

In  principle,  many  methods  for  solving  the  ordinary  differential 
equations  such  as  the  Runge-Kutta  and  predictor  - corrector  methods 
[18]  can  be  applied  to  solve  equation  (9)  with  the  prescribed  boundary 
conditions.  Here,  the  mixed  method  techniques  such  as  the  transfer 
matrix  method  for  beams  are  used.  It  would  appear  reasonable  that 
numerical  integration  could  be  employed  for  shells.  Note  that  the  form 
of  equation  (9)  is  similar  to  the  governing  equation  for  a structural  member 
on  an  elastic  foundation  with  an  equivalent  elastic  foundation  modulus  K 
given  by  equation  (10).  Also,  in  the  above  formulation.  At  must  be  taken 
fairly  small  in  order  to  keep  the  time  discretization  error  reasonably 
small.  However,  a small  value  of  At  corresponds  to  a large  value  of 
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elastic  foundation  modulus  K.  Normally  with  such  a K one  would  expect  to 


I 


encounter  numerical  difficulties  when,  for  example,  the  usual  transfer 
matrix  method  is  applied.  Several  techniques  are  available  to  overcome 
such  numerical  difficulties  [ .2 ].  One  of  the  better  methods  is  to  use 
the  Riccati  transformation  or  the  Riccati  transfer  matrix  method  Cl3- 

Consider  the  formulation  of  the  Riccati  transfer  matrix  method. 
LetCuH  denote  the  NxN  general  transfer  matrix  which  transfers  the 
state  variable  Cw3  from  station  i to  station  i+1  across  segment  i,  at 
time  (n+l)At.  Then 

C wDi+1  - C uDjC W]±  + [ (12) 


where  F contains  the  loading  terms  which  are  evaluated  from  the  general 
loading  function  equation  (11).  Both  W and  F are  Nxl  matrices.  Let  the 
transfer  matrix  be  arranged  and  partitioned  so  that 


Kaite  4M, 


(13) 


where  f contains  the  N/2  state  variables  corresponding  to  the  homogeneous 
left  hand  boundary  conditions  and  e contains  the  respective  complementary 
N/2  state  variables. 

A generalized  Riccati  transformation  at  station  i is  given  by 

Cf}t  - Cs^e^  +C»31  (14) 

which  relates  half  of  the  state  variables  to  the  remaining  state  variables 
at  station  i.  Using  equations  (13)  and  (14),  it  can  be  shown  that  a 
general  recurrence  relationship  could  be  obtained  as  Cl 3 

Cf31+1  - Cs3i+]Ce31+1  +Cp31+1  (15) 


where  [ S]±+1  - [UnS  + U123£U21S  + U^1 


[P] 


i+1  - [U11P  + Vi  - [S]i+l[U21P  + peji 


(16) 


(17) 


Thus,  the  matrices  C S3  andCP3  determine  the  state  variables  f from  e 
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at  any  station.  Another  matrix  [T],  which  transmits  the  state  variables 
e from  station  i+1  to  station  i,  can  also  be  obtained  from  equations 


(13)  and  (14)  as  [1] 

[e]i  " [TWe]i+l  + [Q]i+l 

(18) 

where  [Tl±+1  - [U^S  + U^]’1 

(19) 

[Q]1+l  * -^T^i+i^U21P  + Fe^i 

(20) 

To  start  the  calculation,  the  submatrices  of  the  transfer  matrix 
in  equation  (13)  must  be  determined.  The  matrices  [S],  [P] , [T]  and 
[Q]  are  calculated  for  each  station  while  moving  along  the  member  from 
left  to  right  with  th^  boundary  conditions  [ S ] o * [P]o  = 0.  When  the 
last  station  m is  reached,  equation  (14)  gives 

[f]  » [SI  [el  + [P]  (21) 

m mm  m 

The  N/2  known  state  variables  at  the  right  hand  boundary  are 
substituted  into  the  above  relationship  to  determine  the  remaining 
N/2  unknown  state  variables.  Successive  applications  of  equations 
(18)  at  each  station  allows  the  calculation  of  N/2  state  variables  e 
while  moving  from  right  tc  left  along  the  member.  At  any  station, 
the  remaining  N/2  state  variables  f are  determined  from  equation  (15) . 

This  completes  the  formulation  for  the  continuous  - space  discrete  - 
time  Riccatl  transfer  matrix  method. 

Due  to  the  generality  of  the  transfer  matrix  method,  this  method 
can  be  used  to  solve  structural  members  with  variable  cross  section  and 
arbitrary  boundary  conditions.  In  cases  when  the  transfer  matrix  has  to 
be  evaluated  numerically,  this  method  can  still  be  used  with  the  same 
procedures  as  outlined  above.  For  example,  members  with  continuously 
varying  cross-sectional  properties  and  those  with  more  than  four  equations 
of  motion,  i.e.  N > 4 in  equation  (1),  numerical  integration  may  be 
required  to  evaluate  the  transfer  matrix.  The  finite  difference  discretiza- 
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tion  equations  (2)  and  (3)  are  used  to  transform  the  governing  partial 
differential  equation  into  ordinary  differential  equations,  and  a numerical 
method  such  as  the  Runge-Kutta  or  predictor-corrector  are  applied  to  provide  the 
transfer  matrix  in  equation  (13).  Once  the  transfer  matrix  is  obtained, 
the  Riccati  transfer  matrix  method  can  be  applied. 

Concentrated  occurrences  can  be  treated  the  same  way  as  in  the  Riccati 
transfer  matrix  method  [1]  except  for  a dashpot  or  a concentrated  mass  whose 
behavior  depends  on  the  velocity  or  acceleration  at  that  point.  Consider 
a segment  of  beam  containing  a dashpot  (Fig.  1).  Let  L and  R denote  the 
sections  just  to  the  left  and  right  of  station  i where  the  dashpot  is  placed. 

The  transfer  matrix  across  the  dashpot  can  be  expressed  by  a point  matrix 
at  time  (n+l)At  as 


where  y,  0,  M and  V are  the  displacement,  slope,  moment  and  shear  force  of 

the  beam  respectively.  When  equation  (4)  is  used,  equation  (22)  can  be 
written  in  the  form 


The  same  procedure  can  be  applied  to  the  case  where  the  solution  is 
needed  to  cross  a concentrated  mass.  The  transfer  matrix  across  the  mass 
can  be  expressed  by  the  point  matrix 
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n~  — 
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0 
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Mn+1 

0 

.2 

v 

mh- ±1 

n+1 

L 

3t/ 

(24) 


If  the  approximation  equation  (5)  is  applied, equation  (24)  will  be 
in  the  form 


(25) 


yn+7 

1 0 

0 

0 

yn+l 

"o  “ 

9n+I 

za 

0 1 

0 

0 

9trfl 

+ 

0 

Mn+1 

0 0 
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Mn+i 

0 

Vn+! 

mYo  0 
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1 

Vi 

mQ 

R 

l£pAt2 

L 

n 

The  Riccati  transfer  matrix  method  using  the  point  matrices  (23)  and 
(25)  can  be  used  to  carry  the  solution  across  the  dashpot  and  the 
concentrated  mass.  Such  in-span  indeterminate  conditions  as  moment 
releases  and  rigid  supports  require  special  attention.  The  same  is 
true  for  prescribed  time-varying  displacements.  Consider,  for  example, 
the  beam  in  Fig.  2 which  has  two  successive  prescribed  time  dependent  dis- 
placements 6^  j and  6^  applied  at  position  x » ^ and  a^*  The  unknown 

force  W^_^(t)  is  introduced  at  x • The  point  matrix  corresponding  to 

this  force  is 


(26) 


From  the  condition  ” 6^(tn)  at  t ■ nAt,  x * a^^  and  the  predetermined 
field  transfer  matrix  between  a^  ^ and  a^,  the  unknown  force  ^(t) 
can  be  determined  as 
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<5i(tn)  ■ yi-i(tn)  %y  + ei-i(tn)  uy0  + %m  + vi-i(t  ) uyv 

+ wi-i(tn)  uyv  + uyF 


(27) 


where  U^,  Uy0,  UyM>  UyV  and  Uyf,  are  the  first  row  elements  of  the  field 
transfer  matrix  and  y^_^>  Mi_i  3110  Vi  1 are  c^e  state  variables  at 
x Equation  (27)  then  yields  for  the  unknown  force 


ft  ft  5*? 

yV  yV  yV 


-V 


i-1 


+ ^L(ln)-UYF 


(28) 


U 


yV 


The  point  matrix  to  transfer  across  x = a^  ^ is  found  from  equation  (28) 


as 


1 

0 
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-u  /u  „ 
yy  yv 


0 
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-wv 
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-vv 
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0 
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0 

0 

0 


(6i(tn)-V)/V 


(29) 


This  point  matrix  can  be  used  to  evaluate  the  matrices  C S],  [Pj,  C Tj 
and  CqH  for  the  Riccati  transfer  matrix  method. 

The  investigation  just  completed  is  sufficient  to  find  the  unknown 
force  corresponding  to  an  applied  displacement  using  the  condition  at  the 
next  in-span  or  prescribed  time  dependent  displacement  position.  This 
procedure  can  be  applied  until  the  last  station  is  reached.  However, 
the  N/2  known  state  variables  at  the  right  hand  boundary  are  not  enough 
to  determine  the  remaining  N/2  unknown  state  variables  because  one  of 
the  homogeneous  boundary  conditions  was  used  to  determine  the  point  matrix 
for  the  preceeding  in-span  condition  or  prescribed  time  dependent  dis- 
placement. The  unused  condition  at  the  first,  (counting  along  the  '•member 
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from  left  to  right)  prescribed  time  dependent  displacement  can  re- 
place this  used  boundary  condition  by  being  available  for  fixing  the 
remaining  unknowns.  For  example,  let  i denote  the  station  of  the 
first  prescribed  time  dependent  displacement  and  assume  the  displace- 
ment is  contained  in  the  state  variable  C O*  then,  by  equations  (15)  and 
(18) 


rq.- 

CslTel  +CP1 
mm  m 

(30) 

+ £<0. 
mm  m 

(31) 

Mi’+1 

* ^ T^i+2  ^-e\+2  +^-Q^i+2 

(32) 

Ee]i  ■ 

+^Q\+i 

(33) 

Lf],  - 

CS^Ce^  +CP]i 

(34) 

By  sucessive  substitution  of  equations  (31),  (32),  (33)  and  (34),  the 

following  relationship  between  CfX  and  C e J is  obtained. 

i m 

Cfl  -Cs]  Cel  +CP]  (35) 

l m 

where  [ S]  - [ Cl]±+1  C T]i+2 [T^ 

[P]  - CS],  Cl31+1  Ct11+2 [Q]„  + — 

+ C S]±  C T]1+1  Cq31+2  + c s]j  c q]1+1  + Cp], 

The  prescribed  time  dependent  displacement  of  equation  (35)  and  the 
unused  homogeneous  right  hand  boundary  condition  of  equation  (30)  give 
two  necessary  equations  for  the  two  unknown  right  hand  boundary  conditions. 
Once  the  boundary  conditions  are  evaluated,  the  Riccati  transfer  matrix 
can  be  applied. 

STABILITY  AND  ERROR  ANALYSIS 

how  that  the  formulation  of  the  CSDT  Riccati  transfer  matrix  method 
using  the  general  multistep  method  has  been  presented,  it  is  of  interest 
to  ascertain  the  values  for  p and  the  coefficients  in  equations  (2)  and 
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(3)  to  be  used.  Usually,  it  is  desirable  to  make  p large  so  that  the 
truncation  error  can  be  reduced  C 19 However,  in  structural  dynamics 
problems  wherein  oscillatory  functions  dominate,  it  is  not  clear  that 
the  truncation  error  should  be  used  as  a measure  of  total  accuracy  of 
the  method  E20H.  Some  other  factors  such  as  the  amplitude  decay,  period 
elongation  and  effects  of  spurious  roots  (for  p > 2)  could  also  be 
considered  as  primary  accuracy  measures.  The  stability  and  accuracy 
analysis  of  some  of  the  commonly  used  direct  integration  schemes  in 
linear  structural  dynamics  based  on  the  finite  element  model  of  the 
equations  of  motion  have  been  studied  by  several  authors  [14,21,22,23,24] 
For  our  formulation,  several  techniques  including  the  central 
difference  method,  the  Houbolt  method,  the  Newmark  method  and  the  Wilson 
9 method,  have  been  studied  for  time  discretization  and  it  is  found  that 
Newmark' s generalized  acceleration  method  with  the  following  relations  ■ 
gives  the  most  satisfactory  results 

yn+l  " yn  + At  yn  + (J*  - 8)At2  yn  + 8 ^ yn+l  (36) 

yn+l  " yn  + (1  ' Y>At  yn  + Y At  "y"n+l  (37) 

2 

The  truncation  error  of  this  formulation  is  At  and  no  spurious  roots 
are  involved.  Consider  the  analysis  of  a simple  bar  with  the  governing 
partial  differential  equation 


.2  , .2 

3 y 2 3 y 

3t2  " a 3x? 


(38) 


3t2  “ 3x2 

where  a ■ VEA/m,  E is  the  modulus  of  elasticity,  m and  A are  the  density 
and  cross-sectional  area  of  the  bar, respectively , and  x is  in  the  direc- 
tion of  the  bar  axis. 
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Let  equation  (38)  be  multiplied  by  (4  - 8)At  at  time  nAt  and  by 

2 

BAt  at  time  (n+l)At.  The  following  equation  is  obtained  by  adding  these 
two  expressions  and  using  equations  (36)  and  (37) 

\ 2V  _ « *\  2x 


Ba2At2  3 j^+1  + (h  - 2B  + Y)a2At2  + (%  - y +B)a2At2 


yn+l  - 2yn  + Vl 


(39) 


The  stability  of  equation  (39)  can  be  investigated  using  von  Newmann's 
stability  criteria  C 25].  Assume 

y - T11  eir|x  (40) 

n 

where  T is  a function  of  time  only.  Substitute  equation  (40)  into  equation 
(39)  and  solve  the  resulting  relationship  for  T.  Then,  two  solutions 
are  obtained 


■1,2 


C 2-(^-2B+V)  B2]fcy  D 

2(1+BB2) 


(41) 


where  D = C (V2B+Y)  - 4B(*s-Y+B)IIB  - 4B^ 

2 2 2 2 
and  B^  ■ q a^At 

For  stability,  it  is  necessary  to  have 

D < 0 (42) 

Before  drawing  conclusions  from  equation  (42),  consider  the  error 
associated  with  solution  (40)  by  assuming  that  condition  (42)  is  satisfied. 
Note  that 

, - eil(*±at)  k43) 

is  the  exact  solution  of  the  governing  equation  of  motion  (38).  The 
accuracy  of  the  method  can  be  evaluated  by  comparing  equation  (40)  with 
equation  (43)  and  rewriting  equation  (40)  in  the  form 

(44) 


y - cosn0e  ln(x  ± *a  nAt) 
n 


where  6,  a and  ip  can  be  evaluated  using 

9 _ C (2-XB  )2  + 4(Y+X+B)B2  - (X2-48Y)B4j* 


cos 


2(1+BB2) 


IX,  U 


a = tan 


-1  C 4(Y+X+0)B2  - (X2-40Y)B4j** 


and  X = 20+Y  ; Y = h~Y+& 

From  equations  (43)  and  (44),  we  can  see  that  these  two  solutions  will  be 
identical  to  each  other  if  the  value  of  cos  0 and  ip  are  equal  to  1. 

If  we  set  cos  9=1,  then 

(»S-Y)(0B2+1)B2  = 0 (45) 

This  equation  is  satisfied  if  we  set  y=*j,  i.e.  there  is  no  amplitude 
decay  error  associated  with  the  method  if  Ymh.  Now,  return  to  equation 
(42).  With  Y=*s,  equation  (42)  gives 

D = (1-46)B4  - 4B2 

In  order  to  have  D < 0 for  arbitrary  values  of  B,  we  must  have  0 _>  h » 
i.e.  the  method  will  be  unconditionally  stable  if  we  choose  0 h. 

With  Y”^  and  small  At,  the  parameter  \p  can  be  expressed  as 


\p  = V 4-(  l-40)r|4a4At4 
2-(l-20)n2a2At2 


2 2.  2 

n a At 


V(l-40)n4a4At4 

2-(l-20)n2a2At2 


+ - - (47) 


which  in  no  case  will  be  equal  to  1.  Thus,  the  method  has  a period 
elongation  error  which  is  a function  of  the  time  step  used. 

The  stability  and  error  analysis  for  a simple  Euler-Bernoulli  beam 
can  be  carried  out  in  the  same  fashion  since  the  governing  equation  of 
motion  for  the  beam  has  the  same  form  as  equation  (38)  except  that  there 

are  fourth  order  derivatives  in  x and  EA  is  replaced  by  El.  As  a consequence, 

2 4 

all  n in  the  analysis  should  be  changed  to  n . However,  the  conclusion 

is  the  same,  i.e.  the  method  is  unconditionally  stable  with  0 > h and 
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has  no  amplitude  decay  error  if  y-^j.  Also  there  is  an  error  in  the 
period  elongation  which  is  a function  of  At. 


> 


The  foregoing  stability  and  error  analysis  was  based  on  systems 
with  no  damping.  For  the  case  of  non-zero  damping,  the  stability  and 
error  analysis  would  have  to  include  the  damping  coefficient  as  an  add- 
itional variable.  Usually,  however,  small  values  of  the  damping  coefficient 
do  not  change  the  overall  stability  characteristics  of  an  integration 
scheme  C 3 H . 

It  is  of  interest  to  note  that  if  we  choose  the  parameters  to  be 
8*1  and  y*3/2,  the  Newmark  formulation  equation  (39)  is  identical 
to  the  central  difference  formulation.  We  can  conclude  from  the  analysis 
of  this  section  that  the  central  difference  formulation  used  by  all  of  the 
previous  works  on  continuous-space  discrete-time  method  is  unconditionally 
stable  but  has  amplitude  decay  and  period  elongation  errors.  These  errors 
can  be  investigated  by  setting  8*1  and  y*3/2  in  equation  (44) which  give, 
for  small  At, 


2 2 2 

nQ  -nf|  a At  /2  -,.4. 

cos  0 - e + 0(At  ) 


and  ip  - 1 - n^a2At2/3  + O(At^) 


(48) 


(49) 


Thus,  we  can  see  that  the  central  difference  formulation  of  the  continuous- 
space  di3crete-time  method  has  amplitude  decay  error  approximately  equal 
to  e n n a At  /2  and  the  period  is  increased  in  the  ratio  l/ij>  ■ l+n2a2At2/3. 


NUMERICAL  EXAMPLES 

In  using  the  direct  integration  method,  the  most  difficult  decision 
is  to  select  an  appropriate  time  step  At.  On  one  hand,  the  time  step  must 
be  small  enough  to  obtain  an  accurate  solution.  On  the  other  hand,  the 
time  step  must  not  be  smaller  than  necessary  in  order  to  reduce  the  computa- 
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tion  cost.  Usually  C 3H»  a time  step  smaller  than  0.01T  is  suggested 

for  the  commonly  used  unconditionally  stable  direct  integration  methods 

where  T is  the  period  of  the  response..  Judging  from  the  numerical  results 

of  the  examples  in  this  section,  a larger  time  step  could  be  used  without  losing 

accuracy  for  the  present  method.  However,  in  using  the  Riccati  transfer  matrix 

method  it  is  helpful  if  the  length  of  each  section  is  chosen  to  be  small 

enough  such  that  all  the  terms  for  any  element  in  the  transfer  matrix  are 

the  same  order  of  magnitude  Cl  !•  Hence,  the  section  length  for  the  space 

computations  should  depend  on  the  time  step  chosen  for  the  problem  in 

order  to  satisfy  the  above  criterion. 

It  is  difficult  to  compare  the  efficiency  of  the  present  method  to 

the  other  direct  integration  methods  based  on  finite  element  models. 

The  present  method  is  especially  appropriate  for  structural  members.  The 

cost  of  the  present  method  in  carrying  the  solution  forward  one  time  step 

could  be  estimated  roughly  from  Ref.  1 which  shows  that  the  use  of  the 

Riccati  transfer  matrix  method  to  transfer  from  one  station  to  an  adjacent 
3/2  2 

station  requires  N +4  (N/2)  multiplications,  where  N is  the  order  of 
the  transfer  matrix  for  the  structural  member. 


Extension  bar 

Consider  the  displacement  at  the  tip  of  a cantilevered  elastic  bar 
subjected  to  a suddenly  applied  concentrated  loading  at  the  free  end.  The 
problem  parameters  are  shown  in  Fig.  3. 

The  governing  equation  of  motion  for  this  problem  is  the  wave 
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bar  is  initially  at  rest,  i.e. 


y(„,0)  -0  ; 


t-0 


- 0 


(51) 


To  solve  this  problem  using  the  CSDT  Riccati  transfer  matrix 
method,  replace  the  governing  equation  of  motion  (50)  at  time  (n+l)At  by  the 
following  equation,  in  which  we  have  employed  equation  (36) 


a2y, 


n+1 


dx 

where  R(x) 


Jn+1 

Ba2At2 


' n 

2 2 
(?a^At 


R(x)  ; n-0,1,2,  

y_ 


(52) 


6a2At 


1 A n" 

2 26  ' 1)yn 

a 


Equation  (52)  has  the  same  form  as  a bar  with  an  elastic  foundation  whose 
transfer  matrix  is  C26] 


LPJi+l 


coshai 
EAa  sinh  a l 


sinhai/EA  a 
cosh  a l 


(53) 


2 2 2 

where  a =*  l/f3a  At  and  l is  the  section  length  between  station  i and 
i+1.  With  a fixed  left  end  boundary  condition,  the  submatrices  necessary 
to  perform  the  Riccati  transfer  matrix  method  are 

[f]  - y ; He]  - p 


CUu3  * coshaZ  ; ■ sinha£./EAa 

* EAa  sinh  a£  ; ^C22^  * cosh  <*£ 

l 


(54) 


CFfJ 


R(x)  sinha  (£-x)/adx  ; CFeD  * - R(x)EAcosha(£-x)dx 

Jo 


As  mentioned  earlier,  6 2.  **  should  be  used  in  the  calculation.  In 
order  to  choose  a proper  value  for  6 and  to  see  more  clearly  the  error 
involved  in  the  present  formulation,  suppose  the  bar  is  one  section  long. 
The  displacement  at  the  tip  of  the  bar  for  the  first  time  step  would  be 
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SPoaAtsinhaL 

I?  A t_  ^v,T 


(55) 


The  exact  solution  for  the  problem  is  £27  ] 

* - *8?  * < r (56> 

By  comparing  equations  (55)  and  (56),  we  can  see  that  the  best 
choice  for  tht  value  of  0 is  1 and  in  this  case  the  error  term  due 
to  the  discrete-time  process  is  sinhaL/coshaL.  This  term  will  approach 
unity  as  a + i.e.  At  -*•  0. 

The  exact  and  computed  results  (with  8*1)  for  the  displacement 
at  the  tip  of  the  bar  are  shown  in  Fig.  4.  Several  values  of  8 were 
tried  and  the  solution  does  not  show  much  sensitivity.  Although  the 
present  method  was  proved  to  be  unconditionally  stable,  the  selection 
of  At  used  in  performing  the  integration  determines  the  accuracy  of  the 
method.  Here,  the  accuracy  calculations  were  made  for  different 
values  of  At  and  are  displayed  in  Fig.  5.  These  solutions  indicate  that 
At  = 0.025T  gives  satisfactory  results. 

This  same  problem  is  also  solved  with  damping  coefficient  C * 10  lb-sec/in 
in  order  to  see  the  effect  of  viscous  damping  in  the  application 
of  the  CSDT  Riccati  transfer  matrix  method.  Figure  6 shows  the  results 
of  the  present  method  and  the  exact  solution  found  using  integral 
transforms. 

Euler-Bernoulli  beam  with  simple  supports 

A uniform  Euler-Bernoulli  beam  hinged  at  both  ends  with  a 
bending  moment  applied  at  the  end  x*L  as  shown  in  Fig.  7 is 
analyzed  by  the  present  method.  Numerical  and  exact  C 27D  results  for 
the  displacement  at  the  middle  of  the  beam  are  shown  in  Fig.  8 where  zero 
initial  conditions  are  assumed.  The  numerical  results  match  very  well 
with  the  exact  solution  even  when  the  time  step  0.025T  is  used. 
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The  accuracy  of  the  present  method  as  a function  of  At  for  this 

problem  is  illustrated  numerically  in  Fig.  9. 

In  order  to  see  the  effect  of  amplitude  decay  error  in  using  the 

central  difference  formulation  for  the  present  method,  this  example 

problem  was  also  solved  by  setting  6*1  and  y*3/2.  These  results  are 

2 

shown  in  Fig.  8.  Hartree  C 4U  proposed  the  use  of  the  h extrapolation 
process  to  improve  the  accuracy  of  the  central  difference  formulation  of  the 
CSDT  method.  This  process  was  applied  and  the  improvement  to  the  central 
difference  formulation  is  illustrated  in  Fig.  8. 

Cantilevered  Euler-Bemoulli  beam 
This  example  is  taken  from  Ref.  28  where  the  Newmark  method  was 
applied.  The  problem  considered  was  an  elastic  cantilever  beam  subjected 
to  a suddenly  applied  concentrated  loading  at  the  free  end.  The  parameters 
and  the  finite  element  idealization  of  the  beam  used  in  Ref. 28  are  shown 
in  Fig.  10.  The  numerical  results  for  the  displacement  at  the  tip  of  the 
beam  by  the  present  method  and  the  results  by  the  Newmark  method  C 28]] 
are  shown  in  Fig.  11. 

An  integration  time  step  of  1x10  ^ seconds  was  used  in  Ref.  26. 

The  solution  became  unstable  if  the  time  step  increased  to 

1x10  ^ seconds.  In  the  present  method,  the  beam  was  divided  into  ten 

-4 

sections  and  a time  step  of  1x10  seconds  could  be  satisfactorily  used. 

This  is  about  100  times  larger  than  the  value  employed  in  Ref.  28. 

The  same  type  of  problem  has  also  been  solved  in  Ref.  29  using  the 
Wilson  9 method  with  the  following  parameters  : E ■ 3xl0^psi;  m ■ 0.000733 

4 

lbm/in;  L ■ 120  in,  I ■ 3 in  . Figure  12  shows  the  normalized  displacement 
at  the  tip  of  the  beam. 
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Timoshenko  beam  with  concentrated  loading 
Davids  solved  a cantilever  Timoshenko  beam  by  using  the  so- 

called  direct  analysis  method.  The  parameters  for  the  beam  given  in 

Ref.  30  are:  E = 3xl0^psi;  v ■ 0.3;  y m pg  - 0.3  lb/in^;  A ■ 1 in^; 

2 

A - 0.833in  ; and  IT/I  =*  1,  where  I is  the  transverse  moment  of  inertia 

S L 

and  I is  the  rotary  moment  of  inertia.  Note  that  for  the  Timoshenko 
beam,  the  time  derivative  appears  in  both  the  mass  inertia  term  and  the 
rotary  inertia  term  and  hence,  the  finite  difference  formulation  equations 
(36)  and  (37),  have  to  be  applied  to  both  terms.  The  resulting  equations 
have  the  same  form  as  for  the  static  response  of  a beam  with  equivalent 
extension  and  rotary  foundations  subjected  to  generalized  external  moments 
and  forces.  The  general  forms  of  the  transfer  matrix  for  such  a beam 
can  be  obtained  from  Ref.  26. 

The  numerical  results  from  Ref.  30  and  the  present  method  with  ramped 
shear  and  moment  acting  on  the  free  end  of  the  beam  are  shown  in  Figs.  13 
and  14  respectively.  Only  the  ratio  of  the  transverse  moment  of  inertia  to 
the  rotary  moment  of  inertia  was  given  in  Ref.  30.  No  values  were  indicated. 

A one  inch  square  beam  was  selected  for  our  calculations  and  hence  I ■ I ■ 
0.0833in^.  The  time  step  of  At  ■ 10  ^ seconds  was  chosen.  In  the  direct 
analysis  method,  the  time  step  is  determined  from  At  ■ Ax/C  where  C ■ VEI/pI^ 
is  the  dilatational  wave  velocity  for  the  problem  under  consideration.  The 
section  length  Ax  is  arbitrary  and  is  varied  until  any  reduction  in  this 
quantity  will  not  yield  any  significant  change  in  the  resulting  solution 
of  the  problem.  For  the  parameters  given  in  this  example  problem  C ■ 1.968x10"* 
in/sec.  and  hence  if  30  sections  are  used  as  is  the  case  with  the  present 
method,  the  time  step  will  be  2.54x10  ^ seconds  which  is  much  smaller  than 
the  time  step  used  by  the  present  method. 
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Beams  with  distributed  loading 

The  dynamic  response  of  a cantilevered  rectangular  beam  subjected 
to  a uniform  distributed  ramp  pressure  over  its  length  was  analyzed  by  the 
general  purpose  finite  element  program  MARC  [ 31  ].The  pressure  load 
is  ramped  in  the  first  increment  to  -655.65  psi  and  then  brought  down 
with  constant  slope  to  zero  at  time  .01  seconds.  The  model  with  essential 
parameters  is  shown  in  Fig.  15.  The  beam  was  modeled  with  three  two- 
dimensional,  rectangular  section  beam  column  elements  in  MARC.  Three 
different  methods  of  analysis  were  employed.  They  are  the  Newmark  method, 
the  Houbolt  method  and  the  method  of  modal  superposition.  The  displace- 
ment at  the  free  end  of  the  beam  is  shown  in  Fig.  16.  This  problem  was 
also  solved  by  MARC  using  the  Timoshenko  beam  element  which  allows 
transverse  shear  as  well  as  axial  straining.  The  results  from  Newmark 
and  modal  superposition  methods  are  shown  in  Fig.  17. 

This  problem  is  solved  by  the  present  method  with  a lumped  mass 
model.  The  mass  of  the  beam  is  lumped  at  points  2,3  and  4.  The  point 
matrix  for  a concentrated  mass  with  the  effect  of  rotary  inertia  is 
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where  Fx_  = -pl_(9  +At0  +(Ja-B)At20  )/|3At2 
Ml  T n n n 

and  Fw= -pA(y  +Aty  +(*s-B)At2y  )/BAt2 
V1  n n n 

The  field  matix  for  the  massless  beam  with  the  effect  of  shear  deforma- 
tion is 


1 -i 
0 1 
0 0 
0 0 
0 0 


-£2/2EI 

l/El 

1 

0 

0 


-£3/ 6EI+&/GA 

s 

£2/2EI 


l 

1 

0 


(59) 


where  the  loading  function  F , FQ,  Fw  and  F are 

y y m V 

Fy=p(t)(£4/24EI-£2/2GAg),  F0=  -p(t)23/6EI,  F^=  -p(t)H2/2,  Fv=  -p(t)4 

The  solution  is  carried  out  first  for  the  Euler-Bernoulli  beam  with 
transfer  matrices  (58)  and  (59)  by  setting  1^,  * l/GAg  * 0*  The  result 
is  shown  in  Fig.  16.  The  solution  for  the  Timoshenko  beam  is  also 
calculated  with  the  results  shown  in  Fig.  17. 

Elastic-  lastic  material 

Consider  a cantilevered  bar  with  an  elastic-plastic  material  and  a 

constant  step  velocity  of  100  in/sec  suddenly  imposed  on  the  free  end.  This 
problem  has  been  solved  in  Ref.  28  using  the  Newmark  method  with  a 
finite  element  model  of  the  bar.  The  parameters  used  in  Ref.  28  and  the 
elastic-plastic  model  for  the  materials  together  with  the  finite  element 
model  are  shown  in  Fig.  18. 

The  problem  is  solved  first  for  the  elastic  wave  propagation.  The 
results  for  the  velocity  and  stress  of  the  bar  obtained  from  Ref.  28  and 
the  present  method  are  shown  in  Fig.  19.  The  plastic  wave  propagation 
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problem  is  then  solved  using  the  elastic-plastic  model  shown  in  Fig.  18. 
The  results  from  Ref.  28  and  the  present  method  with  8*1  are  shown  in 
Fig.  20. 

A possible  formulation  for  the  application  of  the  CSDT  Riccati 
transfer  matrix  method  to  non-linear  structural  members  is  proposed 
in  Ref.  32.  Both  material  and  geometric  nonlinearities  are  considered 
and  usually  an  iteration  process  would  be  needed.  However,  for  this 
problem,  advantage  can  be  taken  of  the  fact  that  the  general  stress- 
strain  curve  in  the  elastic-plastic  model  is  composed  of  two  straight 
lines  and,  hence,  the  solution  can  be  carried  out  by  using  the  elastic 
modulus  E to  form  the  necessary  transfer  matrix  if  the  calculated 
stress  on  a particular  section  of  the  bar  is  less  than  the  yield 
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CONCLUSIONS 


The  usual  approach  to  transient  structural’  dynamics  is  based  on 
continuous  time,  discrete  space  ordinary  differential  equations. 

Here  we  propose  a discrete-time,  continuous-space  approach,  with  the 
Riccati  transformation  used  to  stabilize  the  spatial  computations 
for  structural  members. 

The  generalized  Newmark  acceleration  formulation  has  been  chosen 
to  approximate  the  time  derivatives  in  the  governing  equations  of  motion. 

The  resulting  method  is  unconditionally  stable  and  has  no  amplitude 
decay  error  if  the  two  parameters  in  the  formulation  are  chosen  as 
y=4  and  B>Js.  The  method,  however,  has  period  elongation  error  proportional 
to  the  time  step  used  in  the  integration.  Selection  of  the  proper 
integration  time  step  is  always  a problem  in  the  use  of  a direct 
integration  method.  Although  most  of  the  commonly  used  methods  are  un- 
conditionally stable,  there  are  amplitude  decay  and  period  elongation 
errors  associated  with  the  use  of  a large  At.  Usually,  a time  step 
smaller  than  0.01T  is  suggested  for  most  of  the  methods  where  T is  the 
period  of  the  response.  Although  the  same  difficulty  is  experienced  in  using 
the  present  method  for  bar  and  beam  vibration  problems  it  appears  as 
though  a somewhat  larger  time  step  can  be  used  with  sufficient 
accuracy.  This  conclusion,  however,  should  be  scrutinized  after  this 
method  is  applied  to  other  types  of  structures. 

Numerical  examples  included  bar  and  beam  vibration  problems.  Although 
the  geometries  and  loading  conditions  were  simple,  the  proposed  method  ap- 
plies to  more  complicated  situations.  Due  to  the  generality  of  the  transfer 
matrix  method,  the  method  of  this  paper  can  be  applied  to  structural  members 
with  arbitrary  boundary  conditions,  inspan  supports,  and  geometric  and 
material  nonuniformities. 

" - — • «« 
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Fig.  2 Beam  with  prescribed  time  dependent  displacements 
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Displacement  at  x°L  (In) 


x 


U(x.t) 


P(t) 


Po 


a • 0.002  lb-sec^/ln*;  L • lOln;  E • 10?  lb/in2;  A - 1 In2;  P - lO^lb 

o 

Fig.  3 Cantilever  bar  with  a concentrated  loading  at  the  free  end 


Exact  Solution  

Present  Method  S*1  — 

(At  - 0.25T) 


Fig.  4 Results  for  the  vibration  of  the  cantilever  bar  of  Fig.  3 
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Displacement  at  x-L  (in) 


Exact  Solution  __ 

Present  Method  At  ■ 0.025T 

At  - 0.05T  

At  =•  0.1T 


Fig.  5 Effects  of  the  size  of  At  for  the  vibration  of  the  cantilever 
bar  of  Fig.  3 


Exact  Solution 


t (msec) 


Fig.  6 Results  for  the  vibration  of  the  cantilever  bar  of  Fig.  3 
with  viscous  dampling  included 
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Displacement  at  x=I./2  (In) 


Displacement  at  x«L/2  (In) 


0.015 


Exact  Solution 
Present  Method  At«0.025T 


Fig. 


9 Effects  of  the  size  of  At  for  the  vibration  of  the  simply 
supported  beam  in  Fig.  7 


P(t> 


P(t)  ,, 

1 

P 

r 

o 

1"  - 3" 

L, 

0 - 0.000722  lb-sec2/ln4;  E - 3xl0?  lb/in2;  I - 0.832xl0-4  in4;  A - 0.01151  In 


Fig,  10  Cantilever  beam  with  a concentrated  loading  at  the  free  end 


Fig.  11  Results  for  the  vibration  of  the  cantilever  beam  of  Fig.  10 


Modal  Solution 


Fig.  12  Results  for  the  vibration  of  the  cantilever  beam  in  Ref.  29 
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p - 7.675xlO-4  lb-sec2/in4  ; E - 3xl07  psi  ; A - 14.7Zin2  ; v - 0.3 

Fig.  15  Cantilever  beam  with  distributed  loading 


Modal  Solution 

Newmark  Method  _____ 
(6-1/4) 

Houbolt  Method  — — — — 


Fig.  16  Euler-Bemoulli  beam  results  for  the  vibration  of  the 
cantilever  beam  of  Fig.  15 
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Fig.  17  Timoshenko  beam  results  for  the  vibration  of  the  cantilever 
beam  of  Fig.  15 


Exact  Solution 


Newmark  Method 


Fig. 


t(msec) 

19  Elastic  wave  solution  for  the  vibration  of  the  elastic-plastic 
material  bar  of  Fig.  18 


Exact  Solution 


Fig.  20  Plastic  wave  solution  for  the  vibration  of  the  elastic-plastic 
material  bar  of  Fig.  18 
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